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Abstract. The bispectrum vanishes for linear Gaussian fields and is thus a sensitive probe of non- 
linearities and non-Gaussianities in the cosmic density field. Hence, a detection of the bispectrum in 
the halo density field would enable tight constraints on non-Gaussian processes in the early Universe 
and allow inference of the dynamics driving infiation. We present a tree level derivation of the halo 
bispectrum arising from non-linear clustering, non-linear biasing and primordial non-Gaussianity. A 
diagrammatic description is developed to provide an intuitive understanding of the contributing terms 
and their dependence on scale, shape and the non-Gaussianity parameter /nl- We compute the terms 
based on a multivariate bias expansion and the peak-background split method and show that non- 
Gaussian modifications to the bias parameters lead to amplifications of the tree level bispectrum that 
were ignored in previous studies. Our results are in a good agreement with published simulation 
measurements of the halo bispectrum. Finally, we estimate the expected signal to noise on /nl 
and show that the constraint obtainable from the bispectrum analysis significantly exceeds the one 
obtainable from the power spectrum analysis. 
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1 Introduction 

The question whether the inhomogeneities in our Universe have been seeded by a Gaussian initial 
distribution raised a lot of excitement recently (see [1, 2] for reviews). Inflation [3-5] is a theoretical 
paradigm that could generate the initial fluctuations, but has not yet been directly confirmed obser- 
vationally. Standard slow- roll inflation predicts a very low level of non-Gaussianity. However, there 
is no shortage of single and multifield inflationary models with most of them predicting a fluctuation 
distribution distinct from the simple Gaussian case. Thus detection of a non-Gaussian signal would 
provide unprecedented information about the dynamics driving inflation and the interactions of the 
inflaton field [6]. 

The fluctuations in the inflaton field are imprinted in the distribution of photons and matter in 
the Universe. This raises the question, which observable is best suited to detect the tiny deviations 
from the fiducial Gaussian distribution of field amplitudes. Certainly, it is promising to look at 
statistics that would vanish in the Gaussian case, such as the bispectrum of the Cosmic Microwave 
Background (CMB) radiation that for a long time was believed to be the most promising probe 
of primordial non-Gaussianity [7]. In Large Scale Structure (LSS) the detection of primordial non- 
Gaussianity is hampered by the non-Gaussianity produced by late time non-linear clustering, a caveat 
not present in the linear CMB physics. Only in recent years it was realised that equally strong 
constraints can be obtained from the LSS [8]. One of the most promising features of primordial non- 
Gaussianity in the LSS is the scale dependence of the halo bias. This scale dependent bias is most 
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prominent for non-Gaussianities with non-vanishing squeezed hmit of the bispectrum, such as the 
local shape [9, 10] and the shapes with both equilateral and local limit recently found in the Effective 
Theory of Multifield Inflation [11]. It was theoretically predicted for the local type non-Gaussianity 
by [12] and subsequently other derivations were presented by [13, 14]. These were confirmed in 
simulations by [12, 15-20]. First data analysis based on the scale dependent non-Gaussian bias lead 
to remarkably strong constraints on local non-Gaussianity [14]. 

The two point function and its Fourier transform, the power spectrum, are the most impor- 
tant statistics that have been used to analyze LSS surveys so far. Their drawback in terms of the 
detection of non-Gaussianities is that the signatures may be small and difficult to separate from non- 
Gaussianities generated by gravity, whereas the distinct features of alternative inflationary models will 
be imprinted more clearly in higher order statistics such as the three point function or the bispectrum 

{S{k,)S{k2)5{k3)) - {27rf S^^'^ki +k2+ k^)B{k^M. fcs), (l-l) 

which vanish for Gaussian fields. Measuring the higher order statistics is more involved than the 
standard power spectrum analysis because these statistics show both shape and scale dependence, and 
both of these dependencies need to be considered in order to constrain different forms of primordial 
non-Gaussianity. In this paper we try to address whether the bispectrum analysis can improve on the 
power spectrum analysis of LSS surveys in terms of detecting primordial non-Gaussianity of the local 
type. Previous studies of the halo bispectrum in presence of primordial non-Gaussianity [21, 22] are 
based on the standard local bias model. The non-Gaussian correction to the halo bispectrum in these 
approaches arises mainly from loop terms and is thus dependent on the smoothing scale. We will try 
to improve on these calculations and present an independent and smoothing invariant approach to 
obtain the halo bispectrum. 

In [23] it was proven that the optimal estimator of non-Gaussianity is the bispectrum. This raises 
the question how the power spectrum could possibly give tighter constraints than the bispectrum. 
In the analysis of LSS surveys it is extremely difficult to extract information from the very high- 
k non-linear modes. This restricts the analysis of the bispectrum to relatively low-fc modes that 
are sufficiently linear. The situation is changed when one considers biased tracers of the density 
field. These can be related to bispectrum [24], hence they trace non-Gaussian information. The bias 
effectively allows us to extract some of the non-Gaussian information in the high-fc modes from the 
power spectrum analysis. This effectively increases the number of modes in the survey, and allows 
to tighten the limits on the non-Gaussian parameters. Still, one would expect that the bispectrum 
analysis contains further information and the question is how it compares to the power spectrum 
analysis of biased tracers. 

For simplicity, we will focus our attention on the local type of non-Gaussianity, where the poten- 
tial shows a self coupling that is local in real space. The local shape of non-Gaussianity is for instance 
predicted by multi-field infiation [25-27] and in the bouncing cosmology model [28]. Recently [11] 
found new shapes with none vanishing squeezed limit whose LSS phenomenology is yet to be derived. 

This paper breaks down as follows. We first review the basics of non-Gaussianity in Section 
2 and then describe the multivariate biasing scheme and peak-background split approach previously 
introduced by [14, 29] in Section 3. Section 4 introduces the perturbative solutions for the distribution 
of matter and biased tracers as well as our new diagrammatic prescription for the calculation of their 
n-spectra. Section 5 is devoted to the calculation of the halo bispectrum whose constraining power is 
compared to the power spectrum in Section 6. We conclude our findings in Section 7. 

2 Basics 

We consider the local type of non-Gaussianity [9, 10, 30, 31] 

$„G(a;) = v{x) + /nl (<^'(a;) - (^')) + 5NL</'^(a;), (2.1) 

where ip is an auxiliary primordial Gaussian potential.^ Following the peak-background split approach 
[14] we consider the potential as a superposition of small and large scale modes if = ips + separated 

^ The coupling of the potentials in Eq. (2.1) is naturally imposed in the early Universe during Inflation. This 
approach, which is followed by our study, is denoted the CMB convention. However, some authors impose the same 
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by a cut-ofF wavenumber A. Thus from Eq. (2.1) one obtains 



^nG = Vi + hh^i + g^h^i + (1 + S/nlV'i + Sg^hVi )y^s + (/nl + 3.gNLtpi) + 9Nhy^s ' (2-3) 

where all the fields are evaluated at the same spatial position x. Short modes in the above expression 
can be easily identified since only terms containing at least one Gaussian short mode can contribute 
to the short wavelength power. These short wavelength modes dominate the collapse of dark matter 
haloes, whereas the long wavelength modes raise or lower the actual density in large patches of the sky, 
effectively lowering and raising the collapse threshold. In the presence of non-Gaussianity the long 
wavelength modes furthermore affect the variance of the short modes and thus lead to an additional 
dependence of the number density of collapsed objects on the amplitude of the long wavelength modes. 
As we will see, in the presence of local non-Gaussianities, this effect is proportional to the value of 
the long wavelength Newtonian potential, leading to a distinct scale dependent effect. 

The actual effect of the long mode on the variance of the non-Gaussian short modes can be 
estimated as follows. Squaring the short part of the non-Gaussian potential we obtain 

*nG,s = (1 + 4/NL<y5i + 6<?NL<^f + ^J^lV'i + 12/nl5nl</:'? + g.g^LV'f) (2.4) 
of which we can easily compute the expectation over the short modes 

'^ncs = (*'g,s), = (1 + 4./nl</'i + GgNL^i + 4/^L<^f + 12/NLgNL¥'? + ^g^L^f) ctI^,, (2.5) 

where we identified ctq ^ — {^1) s- Here we neglect all correlators of odd number of Lps as well as the 
<7q g = {vt) term. The resulting expression agrees with the expressions previously derived by [12, 14]. 
In contrast to the variance, the three point function or skewness 

{<^lo,s), = (^Inl (v^Dl (1 + 4/nl^i) + 6.gNL {^1)1 (2.6) 

vanishes in the Gaussian case. Similar to the variance, the skewness is rescaled by the long wavelength 
potential. 

Now it remains to connect the non-Gaussian effects on the gravitational potential to the distri- 
bution of matter. In the Newtonian limit, valid well inside the horizon, the Poisson equation relates 
the long wavelength Gaussian potential to the density perturbation, 

(2-7) 

a(k, z) 

where we introduced the auxiliary function^ 

2Pc^D{z)T{k) g{z ^ 0) 
a(k,z)^ —2- — , (2.8) 

which scales as k^/H^ on large scales where the transfer function is unity. As discussed in [32-34], 
on horizon scales unphysical gauge modes and relativistic corrections to the Poisson equation require 
a more careful analysis. Fig. 1 shows the Poisson factor as a function of fc-mode. Note that the 
corrections to the Gaussian spectra are given as powers of /nl/ ct{k). The importance of potential 
and density terms is equal at /c « 2 x 10""* hMpc~^ for /nl = and at fc w 2 x 10"'^ /iMpc~^ for 
/nl = 0(100). Here and in the rest of the paper we use the transfer function for a cosmology with 
iln, = 0.25, f^A = 0.75, as = 0.8, = 1.0. 



equation in the late time evolved Universe (the LSS convention). Therefore one has to be careful when comparing 
quoted constraints on the non-Gaussianity parameters. Namely, the potential evolves as 

ip{x, a) = ^^^ip{x, ao = 1) = g{a)(p{x, ao = 1), (2.2) 
a 

where D{a) is the linear growth factor normalised to unity at ao = 1 and thus g{ao = 1) = 1. In an Einstein-de-Sitter 
Universe the potential is constant in time, whereas it decays as g(a = 0)/g{ao = 1) = 1.34 a: 4/3 in the currently 
favoured ACDM model. 

^Note that we are not writing explicitly the norm of a vector but use the notation k = \k\ 
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Figure 1. Left panel: Poisson factor a(k) relating density and potential via 5{k) — a(fc)$(fc). On scales 
of fc ~ 0.1 /iMpc~^ the potential is smaller than the density by a factor of 10^, whereas they are equal 
on very large scales of k ^ 2 x 10''^ /iMpc'^ The non-Gaussian corrections generally scale as f^h/oik) 
and are thus suppressed on high k's. Note that the Poisson equation in Newtonian gauge receives general 
relativistic corrections as a{k) approaches unity. Right panel: Skewness C3 = a S3 and its first and second 
mass derivatives evaluated for /nl = 1 and (pi = 0. 



3 Bias from the Universal Mass Function 

Galaxies and their host haloes are believed to trace a smoothed version of the underlying distribution 
of dark matter. The relation between the local overdensity in the matter and halo fields is described by 
the bias function. This bias function can be related to the abundance of haloes of mass M described 
by the halo mass function. In this Section we will review the basic Gaussian mass functions and their 
non-Gaussian corrections. Finally, we will present the peak-background split and multivariate biasing 
scheme previously introduced by [14, 29]. 

3.1 Mass Functions 

Haloes are assumed to form at the peaks of the underlying dark matter density field. Numerical 
simulations and analytical calculations indicate that the abundance of collapsed objects can be inferred 
from the distribution of points that exceed the density threshold i5c = 1.686. Following the first studies 
of [35] (hereafter PS) it was found that the mass function, the number density of collapsed objects 
of mass M, reduces to an universal functional form for different redshifts and cosmologies if it is 
expressed in terms of the peak height 



CTnG 



(3.1) 



where a{M) is the variance of the density field smoothed on mass scale M. This definition of i' follows 
e. g. [14], while other authors define v = 5c/ a. All the results presented in this paper can be written 
in terms of either definition by replacing the variables accordingly. Then the number of collapsed 
objects of mass M can be expressed as 

Using a random walk in a Gaussian density field PS derived 
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which is however not in very good agreement with the mass function measured in numerical sim- 
ulations. To improve the agreement, [36] (hereafter ST) proposed a modified version of the Press- 
Schechter mass function 




exp 



' 2 



(3.4) 



where the parameters q — 0.707 and p = 0.3 were obtained from a fit to numerical simulations and A 
is a normalisation factor. 

Performing the random walk using non-Gaussian statistics is more involved, however using an 
Edgeworth expansion of the exponential [37] (hereafter LV) obtained for the mass function 



nLv(M) 
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where the skewness is defined as 
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The above derivation based on the Edgeworth expansion is expected to be satisfactory for low peaks 
only. As LV show in their Appendix B, higher cummulants than ^3 gain importance for high mass 
haloes M > 10^^ h~^MQ. As we will see in the next section, bias parameters are basically derivatives 
of the mass function, and thus the bias parameters derived from the LV mass function should be 
trusted for low and intermediate mass haloes only. The full treatment of the excursion set theory 
using non-Markovian random walks [38] seems to be a promising alternative in the high mass limit. 
In the case of local type non-Gaussianity we obtain for the skewness 



C3 = O'G'S's 



6/1 



NL 



d^g' 



(27r)3 J {2tt) 



:aM{q}aM{q')aM{q + q')Pv<p{q)Pvv{q'), 



(3.7) 



where aM{k) = WM{k)a{k) and Wa/ (fc) is the filter of mass scale M. The LV mass function in 
Eq. (3.5) is not in the form of an universal mass function. We can however rewrite it simply applying 
the chain rule to a S3. 
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where we introduced the auxiliary function 



R{v) = 1 



S^a 
6 



3/2 



- iv 



1/2 



Ids, 



3 dv 



^ ( f 



3/2 _ ^1/2 



(3.8) 



(3.9) 



The presence of the aS-i terms still spoils the universality because this term carries a mass dependence 
via the smoothing scale. As shown in [16] and in Fig. 1, a S3 is only very weakly dependent on the 
smoothing radius, such that we can safely treat it as a constant in the mass function. 

Thus we managed to write the non-Gaussian mass function in the form of an universal mass 
function and can benefit from the known results for universal mass functions. Usually the LV mass 
function is multiplied by a correction factor 7(M) = ngx (-^^ )/'T'PS (-^) to improve the agreement with 
simulations, leading to n^j^iM) = ngT(M)i?(i/). The underlying assumption is that ST corrects PS 
for triaxial collapse and LV corrects PS for non-Gaussianity.'^ Alternative derivations of non-Gaussian 



■^It is not clear how well this statement is theoretically justified, but it is enough for us for obtaining an estimate on 
the non-Gaussian effects on the mass function and the induced biases. 
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mass functions that are valid in various, if not all, regimes have been presented by [38-40]. For our 
study, the important result following from these analyses is that to a good approximation, all of the 
above mass functions can be treated as universal, that is as being function of bd o only, and that the 
inferred values of the biases are not very different. While we are aware of disadvantages of both the 
Gaussian and non-Gaussian mass functions, it goes beyond the scope of this paper to discuss their 
detailed validity or to develop an improved mass function. 

3.2 Multivariate Lagrangian Bias 

The spherical top hat collapse model states that a spherical overdensity collapses to form a gravita- 
tionally bound object once it exceeds a certain overdensity threshold be- Long wavelength density 
perturbations raise or lower the mean density in a patch of the Universe and thus effectively raise or 
lower the collapse threshold b^ ^ b^ — b\. So far, the bias parameters were calculated by expanding 
the local number density in the amplitude of the long wavelength density fluctuation b\ only. As 
we saw above in Eq. (2.5), coupling between long and short modes leads to an enhancement of the 
variance of the short wavelength density perturbations. This enhancement is proportional to powers 
of the long wavelength Gaussian potential and thus suggests to calculate the distribution of collapsed 
objects by expanding the local number density in terms of the long wavelength modes of both density 
and potential [14, 29, 41] 

, \ - dn ^ , . dn , . Id^n^^/ x 2 d'^n ^, , , , 1 d^n ^, . 
n{x) =n + -b,ix) + + -^^b, (.) + ^g^H^W^i^) + 2 O-IO) 

Here we write down explicitly the spatial dependence to highlight the local relation between bias, 
overdensity and potential. The multivariate bias expansion is not stating that density and potential 
are independent parameters, but rather that the different scale dependence of density and potential 
requires one to expand in both of them in order to keep the bias parameters scale independent. 
As we will stress at the end of this subsection, it is important that these fields are restricted to 
long- wavelengths. Notice that the fact that the overdensity depends locally only on b and (f is a 
consequence of the functional dependence of the mass function on b and (p even in the non-Gaussian 
case. Defining bi^(x) — n(x)/h — 1 and identifying the partial derivatives with the bias parameters 
we obtain in Lagrangian space 

bl;ix) = b\M^) + b};M^) + ^fblix) + b^M^M^) + + ■■■, (3.11) 

where bo is the initial Lagrangian overdensity. 

We can now calculate the Lagrangian bias parameters under the assumption that the local 
number density can be expressed with an universal mass function as in Eq. (3.2). 

The presence of a long wavelength mode can be accounted for by replacing be ^ be — bi and 
(Tq fnG in the peak height, such that that the conditional peak height in presence of a long 
wavelength mode can be written as 

f '^c-'^l 

Vf7G(l + 2/nl<P1 + 3gNL<^f + 2/nl<^?) 

and the bias parameters are given by (see Appendix B for details on the calculation) 



(3.12) 
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n dv 
4 d'^n 2 v dn 
n b'i dv^ n b'i dv 
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(3.13) 
(3.14) 
(3.15) 
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=2/nl {5cb\^ - b\„) (3.17) 



6^2 = ^ 2^^^ + - (3.18) 

=4/Sl'5c (fc^o^c - 2feL„) + 6(5c5nl61'o (3-19) 

where all the derivatives are evaluated for i5i = 0, t/Ji = 0. We left the derivatives of the mass function 
unevaluated and obtained expressions that are sufficiently general to enable application to different 
Gaussian and non-Gaussian mass functions."* For instance, the derivatives of the ST mass function 
Eq. (3.4) are given by 



1 9nsT qv — \ p 



nsT dv 2v + iqv)P) 

1 d'^nsT + vpq {qv)'^ — '2qv 



nsT 1/2 (1 + {qv)P) 4i/2 

and for the LV mass function Eq. (3.8) by 

1 dnuY 1 dnsT 1 dR{v) 
n^Y dv nsT dv R dv 

1 d'^nvj 1 d^nsT _ 1 dnsr dR{v) 1 d^R{v) 



(3.20) 
(3.21) 

(3.22) 



ulv dv'^ rigx dv"^ ' ngj^R dv dv ^ R dv^ ^"^ '^'^ 

The amplitude of the halo power spectrum at a given wavelength must not depend on the 
smoothing scale, as it is an observable. Let us imagine to perform the one loop computation of the 
halo power spectrum with two different smoothing scales Ai and A2, with A2 > Ai. Since the final 
answer must not change as we change A, we need to renormalize the bias accordingly. This means 
that the renormalized bias at scale Ai is related to the renormalized bias at scale A2 by a relation of 
the following form [42] 

biS = + (^20^ ^ + &3o) (3-24) 

where tr^ ^ is the real space variance computed including only wavenumbers between Ai and A2. 
Notice that the term proportional to 68/21 arises from having inserted a gravitational interaction 
vertex i<2. This tells us how the renormalized bias parameters change as we change the smoothing 
scale and include loops that renormalize the bias parameters. 

In the PBS method we assume that the long mode is infinitely long. This corresponds to having 
taken the smoothing scale A to be zero (or equivalently infinitely long smoothing length). In this 
regime, loop corrections, that always include only modes that run from zero wavenumber to the 
smoothing scale, are zero by construction, as we took A = 0. This tells us that the PBS method 
provides already the renormalized bias parameters at infinite smoothing length. As we change the 
smoothing length and we make it shorter and shorter, we should include loops and renormalize the 
bias parameters accordingly following formulae similar to the one above. However, the above formula 
is by construction such that the effective bias parameters do not change as we change the smoothing 
scale. Therefore, as we change the smoothing scale, we can simply avoid including the loops, like 
the ones above, that renormalize the bias parameters and use directly the bias parameters that we 
obtain at infinite smoothing length. This will lead to the same answer as if we had included all 
the loop corrections and changed the bias parameters accordingly. Although this statement has not 
been carefully verified from a theoretical point of view, the corrections are anyways quite small on 
large scales and the model seems to be in accord with what is inferred from simulations. In fact 
in [12, 14], a good agreement between peak background split predictions for the non-Gaussian bias 



'Actually one only has to calculate the derivative of uf{u) since the proportionality factors cancel out. 
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boi and simulations has been found. Although there is evidence for a weak dependence of the bias 
parameters on the smoothing scale, it has been shown by [43, 44] that even the second order bias 620 
is quite well approximated by the local bias model in Fourier space. A more precise treatment of this 
point goes beyond the scope of the present paper as it would require calculations at loop level and 
comparison to numerical simulations. 

However, the fact that the peak background split method leads directly to the renormalized 
bias parameters is independent of the assumption of Gaussian initial conditions. Thus it is not 
unreasonable to assume that this fact also extends to the non-Gaussian bias parameters used in our 
study. 

The derivation presented so far assumes an universal mass function, i. e. that all the dependence 
on the long modes is implicitly encoded in the peak height v. However, intrinsically non-Gaussian 
properties of the distribution enter in the non-Gaussian mass function, for example S^a in the case of 
the LV mass function. Following an argument similar to the one we used to derive the ip\ dependence 
of the variance in Eq. (2.5), one can show based on Eq. (2.6) that the three point function in the 
presence of a long fluctuation gets rescaled as {S%j) — (1 + 4/NL<y3i). This dependence on ipi can 
not be encoded in v and thus an additional explicit derivative with respect to the long wavelength 
potential arises, which leads to the following corrections to the bias parameters, for example for the 
LV mass function: 



A601 
Abo2 



2 /nlC; 



3 R 



^ (i/3/2 _ 3^1/2) (3.25) 
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(z/'V2 _ 3^1/2^ + 3 ^^1/2 _ ^_l/2^^ 
(^3/2 _ 3^1/2^ ^ 3 ^^1/2 _ ^_i/2^ 



(3.26) 
(3.27) 



These corrections are generally of the same order as the bias corrections arising from the non-Gaussian 
LV mass function in Eq. (3.22) and (3.23). However, for high v the latter dominate. For realistic 
values of /nl all the bias corrections arising from the non-Gaussian mass function are on the percent 
level. The mass function itself can be trusted at the 10% level only and thus these corrections can be 
safely neglected. 

3.3 Transformation to Eulerian Space 

Observations are performed in the late time evolved Eulerian density field. It thus remains to translate 
the above result to Eulerian space. The halo density fields in Eulerian and Lagrangian space are related 
by [45] 

{l + 5l) = {l + 6){l + 5l). (3.28) 
Finally, one wants to write down the Eulerian analogue of Eq. (3.11) 

5f{x) = b^Six) + b^M^) + ^-^S'ix) + bf,S{xMx) + ^-^v'ix) + ... (3.29) 

In what follows we will absorb the prefactors in the bias parameters, thus 620/2! — 620 and 602/2! — > 
602 • The linearly evolved Lagrangian overdensity can be expanded in powers of the final Eulerian 
overdensity 

(^o = XI = ai*^ + «2'5^ + asS^ + ■■■ (3.30) 

i 

where the expansion parameters are given by the spherical collapse dynamics ai = 1, 02 = 
— 17/21. 03 — 341/567. We calculated the corresponding Eulerian bias parameters using the re- 
lations in [29, 45, 46] 



^10 =1 + 610 (3.31) 

^^0 



6^0 ==2(ai+ 02)610 + a?6^o (3.32) 
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Figure 2. Eulerian multivariate bias parameters for /nl = 100, (?nl = 0. Left panel: Bias parameters for 
S. Right panel: Bias parameters for (p. For (;nl = one has 602 oc /nl) '^oi oc /nl, &11 oc /nl and thus the 
rescaled bias parameters are independent of /nl • bo2 shown by the dashed line shows a pronounced minimum 
for M ~ 1 X 10^* /i~^Mq that is multiplied by /^l and can thus lead to a large contribution. The red lines 
are derived from the ST mass function, whereas the blue thin lines are derived from the LV mass function 
including the explicit Lpi correction of Eqs. (3.25), (3.26) and (3.27). 



^01 =^01 (3.33) 

^02 =bo2 (3.34) 

=^ii«i + &01 (3.35) 

Calculations of the evolved density field are most conveniently done in Fourier space. Thus we 
translate the above expansion Eq. (3.29) to fc-space. 

<5h(fc) = bfoSik) + bfMk) + bi, [5 * S]{k) + bf,[S * if]{k) + * ^](fe), (3.36) 

where we introduced the notation [5 * 8\{k) as a shorthand for the convolution integral. From now on 
we will omit the superscript E and all bias factors should be understood as the Eulerian ones. 

In Fig. 2 we show the mass dependence of the bias parameters arising from the Gaussian ST and 
the non-Gaussian LV mass function. The bias factors of the potential and the density-potential cross 
term are scaled by their /nl dependence. From this figure it is obvious that the second order biases 
can be negative over large parts of the interesting mass range. This can lead to subtle cancellations 
between terms in the resulting galaxy or halo n-point functions. Furthermore, the /^l dependence 
and the pronounced minimum in 6o2(Af) around M — 10^^ H^^Mq can boost the corresponding terms 
on large scales. As noted already by [15] and apparent in Fig. 2, the non-Gaussian mass function 
Eq. (3.8) leads to an additional scale independent offset in the bias. From Fig. 2 we see that the 
difference in the bias parameters is most apparent for high mass haloes, whereas there is only a small 
difference for low mass haloes. The differences between the bias from the two mass functions are very 
small and negligible for our purposes. The same is approximately true even for the other non-Gaussian 
mass functions. Thus we use the Gaussian ST mass function for our numerical predictions. 



4 Perturbation Theory including non-Gaussianity 
4.1 Matter Density Field 

Perturbation theory (FT) aims to solve the cosmological fluid equations using an expansion about 
the linear overdensity 6^\k) [47] 

<5,„(fe) ^ 4i)(fc) + S^^ik) + 5(l\k) + .... (4.1) 
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In our approach, these modes should be thought of as having wavelength longer that the cutoff 
that separates short and long modes. Let us define a linearly evolved primordial density field 
{i.e. evolved without taking into account the gravitational evolution), to be (5m, p- At first order, 
a(fc, z)ip{k) = (5m,p(fc, z), i. e. the primordial Gaussian potential (p is of the same order as the primor- 
dial linear overdensity (5m,p- In the following we refer to these two quantities as first order quantities 
and count the order of terms by counting the powers of first order terms. 

The non-Gaussian self-coupling of the potential introduces non-linearities in the evolved primor- 
dial density field, whereas in the Gaussian case we would have (Sni,p = <5m,p. Transforming Eq. (2.1) 
to Fourier space and applying Eq. (2.7) we obtain for the linearly evolved primordial matter density 
up to third order 

5m,p(fc,2) =«(fc,z)$(fc) (4.2) 

f 

=a{k,z)ip{k) + a{k,z)ft^L J j^-^f{q)(p{k-q) 

+a(fc,z)5NL J J j^ifiqMq'Mk- q - q'), (4.3) 

=4i)p {k,z) + /nl4S ik,z)+ ^nlC'I (fc, ^) ■ (4.4) 

The resulting field is then subject to late-time non-linear gravitational clustering, which introduces 
further couplings. To take this effect into account we use the evolved primordial distribution as the 
source for the late time evolution and insert Eq. (4.4) into the known PT expressions for Siri\ (5m ^ 
and (5m^ (see appendix A) 

s[^lG{k,z)^S^J,l{k,z), (4.5) 



'„G(fe- ^) = / (03^m,U9)^i',Ufe - l)F2{q, k q) 

+/NL<5(?.)p(fc,^), (4.6) 
CG(fe-^)= / (0 / (^4\U9'^)'5ilU9''^)'5ilUfc-9-9''-^)^3(q,g',fc-9) 

+ 2/NL j (03'^m.U9' - 9' ^)^2(g, k-q) 

where F2{ki, ^2) and ^3(^1, k2, fca) are the standard second and third order mode coupling kernels. 
4.2 Halo Density Field 

The halo density including higher order corrections from biasing, non-Gaussianity and non-linear 
clustering can be derived using Eq. (4.1) in Eq. (3.36) 

4(fc) =610 {slZoik) + SlZoik)) + boMk) 

+^20 [S^nG * Cnc] i^) + bn [S^^Ig * ^] i^) + &02 * (k) (4.8) 

=bwS^^l{k) + boi^ik) 

/d'^q f d'^q 

j^(p{q)ip{k-q)+b2o J ^^;^(5Wp(q)(5Wp(fc - q) 

+ho J ^S^J,l{q)Si^lik-q)F2{q,k-q) + aik)h^b,o J ^ifiqMk - q) 

+bn J ^(5«p(g)^(fc-q). (4.9) 
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Note that at lowest order we recover the weh known resuh of [12] 



At this order the potential can be replaced by the density due to the linearity of the Poisson equation 
in fc-space. Thus one finally obtains an expression that is proportional to the density, but non-local. 
The next to leading order, however, is not proportional to 6^{k) since the density and the potential 
are convolved with each other. It is due to this effect that the bias expansion should be performed 
both in S and ip. 

We would like to stress an important subtlety concerning the usage of the bias parameters 
obtained from the peak background split in perturbation theory. Depending on the precision required 
for the calculation, it is possible that higher order corrections from perturbation theory need to be 
implemented. In the diagrammatic description, which is explained in the next subsection, there 
are loop diagrams (convolutions of the fields with some kernel) involving non-linear bias vertices. 
Unfortunately, some of these diagrams are highly dependent on the cutoff A {i.e. the smoothing scale). 
Clearly, physical quantities should not depend on the smoothing scale. This spurious dependence can 
be removed by defining effective or, more precisely, renormalized bias parameters that take into 
account both the tree level and the loop contributions and that are directly connected to observable 
quantities. For example, it is possible to show that loop corrections to the power spectrum originating 
for example from 620 ^nd &30 terms effectively renormalize &10 and &01 [42, 48].^ 

Finally, it should be stressed that even the remaining perturbation theory should be performed 
in a way that ensures that the remaining contributions are independent of the smoothing length 
A~^. This suggests to use a carefully defined perturbation theory, such as for example 'renormal- 
ized perturbation theory' [49], or the recently proposed effective fluid description of cosmological 
perturbations [50]. 

4.3 Diagrammatic Representation 

In the previous section we derived the perturbative expressions for the matter and halo density fields. 
As a consequence of the stochastic nature of cosmological fluctuations, there is no hope to directly 
predict the observed distribution of galaxies and matter in the Universe. Rather, we need to calculate 
expectation values of products of the fields and compare them to the corresponding statistics as 
measured in the sky. The calculation of these statistics turns out to be an involved combinatorial 
task if one goes beyond second order in the fields. 

To facilitate these calculations, we present a diagrammatic representation of the mode coupling 
terms that arise from biasing, non-linear clustering and non-Gaussianity. Similar diagrammatic rep- 
resentations of perturbation theory have been used in the literature [47, 51-53] but we are not aware 
of an intuitive inclusion of all the three effects into one prescription. These Feynman diagrams show 
intuitively which coupling terms arise and can be translated into the corresponding equations by 
straightforward application of a set of Feynman rules. 

Let us start representing the fields, as the basic ingredients of the theory. What we want to 
calculate in the end are correlators of halo or matter density fields, thus we need symbols for the 
outer points, namely Sh and Sm- The latter two are represented by the half filled and filled circles 

^In fact, it is easy to estimate the one loop contribution due to 620 combined with an /jml vertex and an F2 
vertex (see next section for these definitions). This diagram induces an effective bias boi numerically equivalent to 
the one obtained from the peak-background split if one considers small external k's, a very high mass scale and most 
importantly sets the smoothing scale equal to the mass scale of the halo as done in [21]. This is an effect that arises 
from pushing the smoothing scale to very high fc's, where it is unclear whether perturbative calculations can be trusted. 
These renormalized parameters do not depend on the smoothing scale, while the coupling constants with which we 
perform perturbation theory and the contributions from loop diagrams do depend on the cutoff A. This is a behavior 
familiar from quantum field theory. Actually, as mentioned in the former section, the bias parameters inferred from 
the peak background split method can be interpreted as the renormalized ones. Therefore, loop diagrams causing the 
renormalization have been already accounted for. This guarantees that the tree level calculation presented here provides 
the dominant contribution on very large scales. A more precise understanding of this point lies beyond the scope of the 
present paper, where we concentrate only on the tree-level calculation and do not deal explicitly with loop orders. 
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depicted in Fig. 3. Note that when (5m is used for an outer point it always includes all the possible non- 
Gaussian and non-linear contributions up to the considered order, whereas the density as a source 
field is linear. Next, we consider the primordial potential iy9, represented by an open circle. Even 
if a{k)(f{k) = ^in,p(fe), we introduce symbols for both the density and potential to make the 
behavior of the potential terms more obvious and to make sure that the potential terms arise only 
directly from the initial conditions. However, no difference is made between the evolved primordial and 
the late time non-linear matter density field in terms of the symbols, because they can be distinguished 
from the context. For instance, the coupling vertices for gravity are sourced by evolved primordial 
matter fields defined in Eq. (4.4) and lead to non-linear fields (see the discussion of the vertices 
below for more details). Finally, the initial conditions are known in terms of the power spectra of 
fluctuations. Thus we also introduce the power spectrum symbolized by the half filled big circle, 
where subscripts are used to distinguish the density-density, density-potential and potential-potential 
power spectra. 

0<P ©Sm 0(5h ®P 

Figure 3. Symbols use to represent the fields and power spectra. From left to right, the primordial Gaussian 
potential ip, the matter density field Jm, the galaxy/halo density field (5h, and the power spectrum P{k) arising 
from two linked fields. 




Figure 4. Non-linear gravitational clustering leads to a convolution integral weighted by a Fi kernel, which 
we symbolize by a square shaped vertex. From left to right we show the first, second and third order mode 
coupling contributions to the matter density field. Note that F\ = 1. 

The non- linear clustering vertices (we will refer to them also as gravity vertices) Fi , coupling two 
matter density fields, are represented by the open squares shown in Fig. 4. The input density field on 
the right hand side is an evolved primordial field defined in Eq. (4.4) , which receives higher order 
corrections from non-Gaussianity only. Fig. 5 shows the vertices arising from the multivariate bias 
expansion, where the input can be a matter density field of any order or the primordial potential. The 
identification of the higher order bias terms with vertices is possible, since in fc-space the products 
of fields lead to convolution integrals, which are similar to the non-linear gravity terms, where the 
Fi kernels are replaced by the scale independent bias factors bij. Thus we interpret biasing as an 
unweighted convolution of source modes. As noted above, loop diagrams involving higher order bias 
vertices can thus lead to large or divergent contributions, effectively renormalising lower order bias 
parameters. Similarly, the non-Gaussian terms are effectively coupling potential modes to a higher 
order primordial matter mode 5^}^. This interaction is represented by the diamond shaped open 
vertices depicted in Fig. 6, corresponding to a(fc)(5NLi where Qnl — {gnL: /nLi • ■ ■}• These coupling 
vertices are a specific property of local type non-Gaussianity, basically expanding the bispectrum in 
terms of the primordial potential 

B$$<t,(fci,fc2,fe3) = 2h'LP^^{ki)P^^{k2) + 2 eye. (4.11) 

For other shapes of non-Gaussianity such an expansion might take a different form with a fc-dependent 
kernel that can be easily included. 

Time integration is trivially performed in standard perturbation theory since all the initial fields 
are linearly evolved, considering only the growing mode. Thus the lines correspond to propagators. 
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i. e. linear growth factors. However, these Unear growth factors can be used to transform the initial 
fields to linearly evolved fields such that we can totally ignore the time evolution as long as we use 
linearly evolved primordial matter fields as a source. We only need to ensure causality by following 
the time evolution of the fields. Non-Gaussian coupling happens directly after inflation, then non- 
linear clustering and finally biasing. The lines used to connect the primordial potential with the 
non-Gaussian vertices and the non-Gaussian bias vertices are dashed to highlight the fact that the 
coupling of primordial potentials, and thus the imprint of non-Gaussianity, happens directly after 
inflation. To facilitate the distinction of the density propagators, we use straight and wiggly lines for 
the matter and halo density field, respectively. 

The i-th order contribution to an evolved matter or halo field can be obtained following the time 
evolution step by step starting from the initial conditions and going all the way to the final field 

1. Draw i initial fields (Sm.p or ip. 

2. Draw the non-Gaussian Qnl vertices and connect them to the primordial potential using dashed 
lines. 

3. Draw the gravity vertices and connect them to the non-Gaussian vertices or initial density fields 
by solid lines. 

4. For biased tracers, draw the bias vertices and connect them to either initial density fields, gravity 
vertices, non-Gaussian vertices or primordial potentials. Use a wiggly line to connect them with 
the outer point. 

So far we focused our attention on the fields. To compute the i-th order contribution^ to the 
n-spectra we need to glue n diagrams with i source fields and n outer points in all possible ways and 
then pair the source fields in all possible ways. Two linked source fields lead to a power spectrum 
and a momentum conserving delta function {6{k)S{k')) = (27r)35(^)(fc + k')Pss{k). The potential- 
potential and density-potential power spectra are related to the density-density power spectra by the 
appropriate Poisson factors a{k,z). For the translation of the above diagrams into mathematical 
expressions we assign fc-vectors to each of the outer fields. The different ways of performing this 
assignment are accounted for by the cyclic permutations of the fc-vectors. 




Figure 5. The multivariate bias expansion can be expressed by the triangular vertices shown above. The 
interaction of two fields on a vertex corresponds to a convolution integral over the ingoing fc-vectors without 
any weighting. The vertex is connected to the outer point by the wiggly halo propagator. Ingoing potentials 
are always primordial because the coupling of long and short modes, and thus the enhancement of the short 
wavelength variance, happens in the early Universe. 

4.4 Feynman Rules for the ri-Spectra 

Even though the diagrams can be straightforwardly translated into the corresponding mathematical 
expressions we write down the Feynman rules explicitly for the sake of definiteness. For the calculation 
of the i-th order contribution to the n-spectrum do the following 

®The index i has to be even since the correlator of an odd number of Gaussian fields vanishes due to the Wick-theorem 
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Figure 6. The non-Gaussian vertices correspond to an unweighted convolution integral times one factor of a 
and the corresponding Qnl = {cnl, /nl, • • •}• We introduce the bnl vertex to represent the conversion from 
(fi to (5m,p- Note that the ingoing potentials are always primordial, which is highlighted by the dashed line. 



1. Draw all distinct connected diagrams with n external lines up to the desired order i in ip 



lU, 

iv, 

V, 

vi, 

vii, 
viii, 
ix. 



For each vertex with ingoing momenta Qi and outgoing momenta pj write a delta function 

(E,9.-E,p.) 

Assign a linear power spectrum {27r)^ 5'^-^\q + q')Piiniq) to each of the big circles with 
outgoing momenta q and q' . Divide by a{q) or a'^{q) for Ps^p and P^^p, respectively. 

For the outer fields with momenta fej write a delta function 6^^^ ki) 

For each square shaped vertex Fn with ingoing momenta qi write a mode coupling kernel 
Fr,{qi,...,qn) 

For each triangular shaped vertex write a bias factor bij 

For each diamond shaped vertex Qnl — {gnl, /nl, 5nl, ■ • •} with outgoing momentum q 
write a(q)QNL 

Integrate over all inner momenta J d^qi/ {2Tr)^ 

Multiply with the symmetry factor 

Sum over all distinct labelings of the external lines 



2. Add up the resulting expressions from all diagrams 
4.5 Matter Power Spectrum 

As a first application of our diagrammatic approach we write down the terms contributing to the 
non-Gaussian matter power spectrum and show the corresponding diagrams in Fig. 7 



1 

(k) -Pll(fc)+P22(fc)+Pl3(fc) 

d^q Piq) P{k-q) 
d^q P{q)P{k-q' 



^4a(/c)/NL 
P{k) 



/NL- 



(27r)3 a{q) a{k — q) 
d\ P{q) 



F2{q,k~q] 



(4.12) 



a{k) J (27r)3 a{q) 



F2{q,k - q)a{k - q) 



P{k) f d^q P{q) 



6a(A:)gNL- , 2/ ^ 

a(k) J [2tt)'^ a (9) 

where P13 and P22 are the standard one-loop corrections to the power spectrum (see [47] and Appendix 
A). The subscripts on the brackets in the above equation can be used to identify the corresponding 
terms in Fig. 7. The functional form agrees with previous results as presented in [54], who also 
performed a numerical evaluation which is thus not repeated here. The corrections arising from the 
/nl terms are generally small and most prominent for high k, where the validity of perturbation 
theory has to be doubted. 
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Figure 7. Diagrams contributing to tlie matter power spectrum as calculated in Eq. (4.12). The first line 
shows the standard PT terms, whereas the terms in the second and third line are purely non-Gaussian. 



4.6 Matter Bispectrum 

A non-vanishing matter bispectrum beyond the non-linear gravitational contribution would be a 
direct sign of non-Gaussian initial conditions. From the diagrams depicted in Fig. 8 we can derive 
the following expression for the tree-level matter bispectrum 



,{k^M. fcs) = \2P{ki)P{k2)F2{k^M) + 2 eye. 



2/i 



P (fci)P(fc2)a(fc3) 
a{ki)a{k2) 



2 eye. 



(4.13) 
(4.14) 



Here eye. is used to symbolise that the function arguments in the preceding terms have to be cyclically 
permuted as {{ki, k2, k^), {k2, k^, fci), (fca, ki, k2)}- An evaluation of the latter term and comparison 
to n-body simulations is provided in a recent study by [55] . They find that the inclusion of one-loop 
terms leads to a considerable improvement of the agreement between theory and simulation on scales 
exceeding fc = 0.1 hMpc~^ . 

Estimating the bispectrum of the dark matter distribution is an involved task even for the 
upcoming lensing surveys. Thus we will focus our attention on the bispectra of biased tracers such as 
galaxies in the next section. 




Figure 8. The two diagrams contributing to the tree level matter bispectrum Bmmm in Eq. (4.14). 
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5 Bispectra of Biased Tracers 



In the following section we present the main result of this paper, the derivation of the halo bispectra. 
The bispectrum measures the correlation between three fields and thus halo and matter fields lead 
to four possible combinations. The matter auto-bispectrum was already discussed in the previous 
section such that we can focus our attention in this section to bispectra involving at least one biased 
tracer. 

5.1 Halo Bispectrum 

Summing up the diagrammatic expressions shown in Fig. 9 we can write down the tree-level expression 
for the halo auto-bispectrum 

i3hhh(fei, fc2, fca) =blo (2P{h)P{k2)F2ik^,k2) + 2/NL ^^^'\\^i^'2"i^'^ + 2 eye 



A 



B 



C 



a{ki)a{k2) 

+*;A,(2P(k,)Pfe)(^ + ^)F,fe,*,) 

a[ki)a[k2) \ct[ki) a{k2) J 

+''10^01 2— — — — — F2 fci,fe2 +2/nl 2fi. \ 2/1. \ ^2 eye. 

V a(ki)a{k2) a''{ki)a^[k2) 

V a^{ki)a^{k2) / ^ 

( P(ki)P(k2) f I 1 \ 

+boibiobo2 2 ^ ) ^' — - + — - + 2 eye, 
V a(ki)a{k2) \a[ki) a[k2) J 

+^?o^o.(2^§l|^-^2cyc.) (5.1) 

V a(ki)a{k2) /p 

u r. fP{ki)P{k2) f I 1 A o 



\a{ki)a{k2) \a{ki) a[k2) 



-bM,, [ Pik,)Pik2) + + ^ ) + 2 eye. 



H 



V a(ki)a[k2) 

+boibiob2o (2P{k,)P{k2) f-^ + ) + 2 eye, 

V \a[ki) a[k2) 

+6^0620 (2P(fci)P(fc2) + 2 eyc.)L 



K 



The subscripts on the brackets can be used to identify the terms with the corresponding diagrams 
in Fig. 9. In the diagrams and the above equation there are no explicit ^nl terms, since the 
vertices are not explicitly contributing to the tree level halo bispectrum. However, there is an implicit 
dependence via the second order bias parameter 692 • We will not further examine this dependence since 
it changes only a parameter of the model. Note that for ^nl — the bias factors scale approximately 
as 5oi cx /nl, &02 oc ^'^d 611 oc /nl- Thus the exponent of /nl in the terms is the same as the 
exponent of a{k) in the denominator. Hence, it is the ratio of f-Mh/ctik), where k = minj/ci, ^2, ^3} 
that is dominating the overall amplitude on large scales. The highest order contribution is /^^ in the 
D term. This term is negative and gains importance for extremely small k and large /nl- 

The wave vectors are related by fci -1-/22 + ^3 =0 and thus the configuration is fully determined by 
the magnitude of one vector fci, one angle = fei • fc2/(fcifc2) and the ratio of two vectors X2 = k2/ki. 
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The magnitude of the third vector is then given by 



X3 = ^ - ^l + 2^a;2+x2, (5.2) 



which for the isosceles configuration fci = ^2 simplifies to X3 = -\/2(l + fi). The transfer function 
is unity for large scale modes entering horizon after matter radiation equality and is damped on 
small scales, leading to an asymptotic slope of ~ —1.75. Starting from a primordial power spectrum 
Po(fc) = ^fc" with n « 1 we see that the primordial Gaussian potential is scale invariant with 
Pipip cx k^^. The matter power spectrum in contrast is given by Pgs = T^{k)PQ{k) and thus it scales 
as k^ at low fc's and approximately as k~^-^ at high fc's. For the isosceles configuration and for low 
k = ki = k2 the dominant contribution to terms including second order bias (terms D-L of the above 
equation) scales as 

a'(fci)aJ(fc3) A:2«+2j2.2j 3 . 

while the dominant contribution to the A, B, C terms scales as 



Pikl)P{k3) ^ k^X3 ^ fc4-2»-2j^l-2, 



Q;*-l(fci)Q;^(fc3) /fc2i+2i- 



■X 



2j ~ -^3 
3 



(5.4) 



The latter equation considers only the /nl contribution, dominating at small X3's. In Table 1 we 
quote the exponents i,j and the corresponding power of /nl- The combination 2(i + j — 1) is the 
exponent of the dominating short mode in the squeezed limit and can go up to k^. The estimation of 
the importance of the terms is further complicated by the different bias prefactors. We thus evaluate 
the expression numerically and discuss the results in Sec. 5.3 below. 

Recent studies of the tree-level bispectrum [56] are based on the univariate bias parameters biQ 
and 620 only (univariate or Gaussian biasing). At tree level they thus consider only a subset of the 
above terms (see their Eq. (18)) leading to 

Bhhh(fei, fe2, fc3) - [BF,{ki,k2, ks) + fc2, ^3)] + &?o^20 [2P(fci)P(fc2) + 2 eye] , (5.5) 

where Bp^ and -B/^l are implicitly defined in Eq. (4.14). This approach neglects the influence of non- 
Gaussianity on the bias parameters, whereas it has been explicitly shown in simulations [15, 18, 29] 
that local non-Gaussianity introduces a scale dependent bias. 
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Figure 9. Tree level bispectrum diagrams considered for Eq. (5.1). The scaling of the different terms with 
/nl as well as their k dependence are described in Table 1. 
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Table 1. Order of the terms contributing to the bispectrum Bhhh in Eq. (5.1) as defined via Eq. (5.3). From 
left to right we quote the exponents of the fc-vectors, the exponent of the long fc-mode, the exponent of the 
long-short ratio and the exponent of /nl- More than one line per diagram can arise if there are different 
possibilities for combining the components of the diagram. 



5.2 Cross-Bispectra 

The measurement of the halo auto-bispeetrum is limited by the finite number of haloes and shotnoise. 
Thus it is interesting to eonsider also the cross spectra between matter and haloes. Weak gravi- 
tational lensing in cross-correlation with galaxies could be a possible observational probe providing 
measurements of these effects. Without showing the diagrams we write down the expressions for the 
halo-halo-matter bispectrum: 

i?hhm(fci, fc2, fca) =bio (GF2{k,,k2)Pih)Pik2) + 6/NL«(fc3)^^^^ + 2 cyc. 

\ a[ki) a{k2) 

1 1 



+601610 4F2(fci,fe2)P(fcl)P(fc2) 



a{ki) a{k2) 



+^hLaik3)P{ki)P{k2) ( \ + \ ) + 2 cyc. 

\a-^{ki)a[k2) a(ki)a^[k2) J 

\ a{ki)a[k2) a^[ki)a^[k2) 



+610620 (^4P(A;i)F(fc2) + 2 cyc. j (5.6) 
+620601 (2P{ki)P{k2) f — ^ + + 2 cyc. 

+611610 (2P{ki)P{k2) i^r^ + + 2 cyc. 
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a'^{ki)a'^{k2) 



1 




a{ki) 


«(fc2)) 
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«(fc2)) 


a{ki) 
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2 




a{ki)a 


+ 2 eye. 


) 



-611601 ( P{ki)P{k2) y-z^T7^^ + „ X + T^TTTT ) + 2 cyc. 

P{ki)P{k2) 
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+602601 2P{ki)P{k2) . . + . . + 2 eye. 

The eorresponding predietion of a purely univariate bias model is given by 

Bhh^{ki,k2, fea) = 36^0 [BF,{ki,k2,k3) + Jfei, fcs, feg)] + 2610620 [2P(fci)P(fc2) + 2 eye.] . (5.7) 
Finally, one ean also eorrelate two matter and one halo density field: 

-Bhmm(A;i, fc2, fca) =bio (m{kuk2)P{ki)P{k2) + 6/NLo(fc3) '^if'j'^.^f"^ + 2 eye. 

\ a{ki)a[k2) 

1 1 



+601 2F2(fci,fc2)P(fcl)P(fc2) 



a(ki) a{k2 



+2f.Mh)Pik.)Pik.) (^^..(.^ + ,^(,^ ) + 2 eye. 



+62o(^2P(fci)P(fc2) + 2eyc.j (5.. 
+6n(p(fcOP(fc2) (^ + ^)+2eye. 



P{ki)P{k 



2 



2 eye. 



a(^i)a(fc2) 

where the eorresponding predietion of the univariate bias model reads as 

Shmm(fci, fc2, fca) = 36io [Bp, (fci, fc2, ^3) + P/«Jfei, ^2, fc,3)] + 620 [2P(fci)P(fc2) + 2 cye.] . (5.9) 

As for the halo auto-bispeetrum, we show the results of a numerieal evaluation of the above results 
in the next subseetion. 

5.3 Discussion of the Results 

In this seetion we foeus on the evaluation and diseussion of the bispeetrum expressions derived in Sec- 
tion 5 above. For the evaluation no smoothing is required sinee at tree level there are no integrations. 
The effeet of smoothing at tree level breaks down to a multiplication with a smoothing function and 
thus suppresses the result on small scales keeping the large scales unaffected. The bias parameters 
are evaluated for a peak height oi v = A corresponding to haloes of M w 1 x 10^'' h~^MQ and our 
fiducial cosmology introduced in Section 2. 

In Fig. 10 we evaluate the terms in Eq. (5.1) separately to asses their importance. To make the 
comparison of the terms easier we plot the reduced bispeetrum 

^^'^''^''^^ ^ pik,)Pik2) + Pik2)Pik,) + Pik,)Pik,y (^-^o) 

• We see that diagram A's contribution is constant at high-Zc's and grows at low-Zc's. This is due 
to the fact that at high-fc's the signal from the subdiagram Ap^ grows (a fact that by the way 
Q is designed means that Q is constant) , while on large scales the signal from the A f^^^ grows. 
The two subdiagrams have a different scale dependence. 

• It is easy to see that diagrams B are obtained by substituting (5-bias with the one in ip. Sinee 
if is suppressed with respect to S for high-fc's, we have that the signal in B is peaked at smaller 
fc's. The value of /nl used for this plot is rather large, thus the B contribution is larger than 
the one of A at small fc's. 

• In order to pass to the C diagrams, we exchange one 6 leg with a ip leg and exchange a 6io factor 
with a 6oi. Since in the B diagram one of the legs is already associated to the low fc mode, 
we need to associate the remaining leg to a high k mode. This means that the C diagram scales 



- 20 - 




k[hMpc"i] k[hMpc"^] 



Figure 10. Scale dependence of the terms contributing to tlie lialo bispectrum Bhhh in Eq. (5.1) for /nl = 100 
and 2/ = 4 in an isosceles configuration k = ki = k2, ^l = —0.99. The lines are labeled according to the 
corresponding diagrams in Fig. 9 and Table 1. The thin red dashed line in the left panel shows the bispectrum 
in the Gaussian case, which vanishes on large scales due to the shape dependence of the F2 kernel. We show 
the reduced bispectrum defined in Eq. (5.10). The amplitude of the terms is determined both by the bias 
prefactors and the amplitude of the fc-dependent factors. We see that the B and H terms are dominating on 
intermediate scales A; ~ 1 x 10~^ ftMpc"^, whereas D,E and G take over on smaller fe's. 

with respect to the B diagram by a factor of 6oi/(&ioQ^(fchigii)) that, for the values of /nl and 
fchigh, that we are using is about a factor of 10 at the smallest fc's we plot. Notice that since the 
terms here have additional factors of 1/fc^, it raises more steeply at low /c's. Even more than 
for the B diagram, this term becomes irrelevant at high fc's because of the many factors that 
it involves. 

• The remaining diagrams D — L are constructed using the non-linear biases in all possible con- 
tractions. This means that the steepest at low fc's will be the D diagram that involves 602 and 
b'^i, while the most important at high fc's will be the L diagram, that does not involve any ip. 
By construction the L diagram is scale independent when plotted in terms of its contribution 
to the reduced bispectrum Q. 

In Fig. 11, we show how the various diagrams scale with the shape parameter xa in the isosceles 
configuration. For to overall size we consider fc = 0.1 /iMpc~^ and fc = 0.01 /iMpc"^. The scaling 
is also given in Table 1 and can be easily reconstructed from Eqs. (5.3) and (5.4). There is a clear 
correction at small fc's from the Lp terms. Since the signal is not scale invariant, as we move to higher 
fc's, the contribution from terms involving and its non-linearities, such as the Ap^ diagram, gain 
importance. 

The left panel of Fig. 12 shows the dependence of the bispectrum amplitude on the variation 
of the non-Gaussianity parameter /nl- This plot shows that higher than linear powers in /nl have 
to be considered to describe the full bispectrum amplitude. The plot also shows that the sensitivity 
to /nl increases as one considers more squeezed configurations. While evaluated for a different 
halo sample and redshift, the plot shows qualitative agreement with the simulation measurements of 
[57]. In the right panel of Fig. 12 we plot the angular dependence of the bispectrum for an almost 
isosceles configuration fci = fc2/1.2 since the exact isosceles configuration is divergent for 9 tt. Here 
6 = arccos (/.t) refers to the angle between fei and fc2-^ Besides an overall enhancement the most 
remarkable feature is the upturn of the multivariate bias prediction in the squeezed limit. 

'^The divergence of the isosceles case means nothing but the fact that as 9 — > tt one of the sides of the triangle goes 
to zero, and so the bispectrum amplitude goes formally to infinity. 
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Figure 11. Dependence of the terms contributing to tlie lialo bispectrum Bhhh in Eq. (5.1) on the long-short 
axis ratio 13 for an isosceles configuration with k = k\ = k2- The overall scale of the triangle is fc = 0.1 /iMpc~^ 
for the upper panels and 0.01 /iMpc~^ for the lower panels. This plot clearly shows the effect of the (p-hias 
corrections in the squeezed limit 2:3 0. As in the above plot the thin red dashed line in the left panels shows 
the contribution. 



Fig. 13 shows the halo auto- and cross-bispectra. We show both their scale dependence as well as 
the ratio to the Gaussian expression for a squeezed isosceles configuration fei = ^2, = —0.99, ~ 1/7 
and compare to the predictions of the univariate bias model. These figures show that the non-Gaussian 
bispectrum asymptotes to the Gaussian one for high fc's (small scales). On larger scales both the 
univariate bias model as well as our prediction are increased with respect to the Gaussian case due to 
the coupling between long and short modes. Our expression for Bhhh predicts an enhancement by a 
factor of 2 with respect to the existing univariate biasing calculations on scales of fc « 0.03 hMpc~^. 
The enhancement is less pronounced for the mixed halo-matter bispectra. From the scaling with k 
and 2:3 in Eq. (5.3) and Table 1 it is clear that the bispectrum amplitude is largest for squeezed 
configurations /i — —1, 9 tt and low overall scale A: — 0. However, the shape and scale of the 
triangle are limited by the fundamental mode and the sampling variance in exactly this limit. We 
will quantify the signal-to-noise ratio in the next section. 
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Figure 12. Left panel: Dependence of the amplitude of tiie lialo bispectrum -Bhhh on tlie non-Gaussianity 
parameter /nl for (;nl ~ and fci = ^2 = 0.042 /iMpc~^. We sliow the bispectrum amplitude for three 
different values of the short-long ratio defined in Eq. (5.2). This plot agrees qualitatively with Fig. 4 
in [57]. It is obvious from this plot that the bispectrum decreases with respect to the Gaussian non-linear 
clustering case for weakly negative /nl. Right panel: Angular dependence of the non-Gaussian and Gaussian 
reduced bispectra for /nl = 100 as a function of ^ = arccos (fei, k2) — arccos (fi) for a halo sample with i/ — 4 
and a nearly isosceles configuration with ki = ^2/1. 2 = 0.03 /iMpc~^. We show our result (5.1) (red dashed) 
and the univariate biasing result (5.5) (blue dash-dotted). Besides the overall enhancement of the multivariate 
bias result with respect to the univariate bias result there is a remarkable upturn for 9 n. 



5.4 Comparison to Simulations 

While a full simulation based analysis of the halo bispectrum goes beyond the scope of this paper, we 
can nevertheless compare to published simulation measurements. The only published simulation mea- 
surement of the halo bispectrum we are aware of, was performed on a L = 2000 h^^Mpc cosmological 
simulation by [57]. They consider a cosmology with erg — 0.816, — 0.28, Ha — 0.72 and present 
their results for z = 0.5. In Fig. 14 we show the comparison between our theoretical predictions and 
their simulation measurement for a halo sample with Af > 4.6 x lO^'^ H^^Mq. The bias value inferred 
from the simulations 6io,sim — 2.9 [58] points towards a halo sample with M > 7 x 10^^ h~^MQ. This 
deviation to the low mass cutoff quoted by the authors might be due to an incomplete sampling of 
haloes around the cutoff. We adopt M > 7 x 10^^ H^^Mq as the lower boundary of the halo sample 
and calculate the average bias parameters using the LV mass function. Consistent with [57] we rescale 
&01 as &01 ~^ 0.75 &01: a- modification motivated by ellipsoidal collapse (although this seems to apply 
only to Friends-of- Friends identified halos and not to spherical overdensity halos, [2]). Fig. 14 shows 
that the scale and shape dependence is quite well described by the theory, especially the upturn of the 
bispectrum for squeezed configurations. The terms dominating the non-Gaussian signal for « 0.1 
in this comparison are A,B,H and I in Eq. (5.1). Concentrating on the dominating terms for the 
squeezed isosceles configuration ki = k2 ^ k, the bispectrum reads as 



-Bhhh(fei, fe2, fes) = (46?o/nl + 2610^11) I , ^ + (4/NLfo?ofooi + 26io6iiboi) 



ayx^k) a'^[xj,k) 
+ (Af }?h +Ah h h ^mPix^J^..2, P{k)P{x,k) li-^xl 

+ (4/NL&10&01 +4610611601) «(^3fc)«(fc) + ^10^01 ^(^ ^ (5.11) 



- 23 - 



10 



10 



10" 



- 10" 
10' 
10^ 
10^ 

10'^ 

10' 

10" 

\ 10' 
2 

CO 

^ 10" 
CQ 

10' 
lO'' 
10= 




10 



10' 

k [h Mpc ' 



10 



-<NL=0 

f|^i_=1 00 univariate bias 



10 



10"' 

k [h Mpc"''] 



10 



'nl=0 

fiji =1 00 our prediction 



f|gi^=100 univariate bias 




20 



f|^i^=100 univariate bias 



10^ 


10^ 

k [h Mpc"" 


10° 

] 








f|^i^=100 our prediction 


20 


"\ 

\ 




f|^i_=1 00 univariate bias - 




\ 

\ 






10 


\ 






8 


\ 








\ 






6 


\ 








\ 






4 


\ ^ 








N ^ 








\ 








^. X 






2 


^. N 
N \ 








S \ 










N 






"S. 














10-= 


10" 


10° 



l< [h IVIpc'] 



20 



10 
8 

o 

m 6 

o 



. f|vjL=100 our prediction 
f|^i_=100 univariate bias 



\ 

\ 

' \ 
\ \ 
\ \ 



10 

1^ [h Mpc"'] 



10 



Figure 13. Upper Left Panel: Halo auto bispectrum Bhhh according to Eq. (5.1) for /nl ~ 100 and /i — —0.99 
in the isosceles configuration k\ = k2 = k. We show the Gaussian bispectrum (black solid line), the univariate 
bias prediction of Eq. (5.5) (blue dash dotted) and our predictions (red dashed). Upper Right panel: Ratio 
of the non-Gaussian bispectra and the Gaussian bispectrum shown in the left panel. Middle Panels: Same as 
above, but for the halo-matter cross bispectrum Bhhm according to Eq. (5.6). Lower Panels: Same as above, 
but for the halo-matter cross bispectrum Bhmm according to Eq. (5.8). 
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Figure 14. Halo bispectrum in the isosceles configuration fci = for z = 0.5 and M > 4.6 x 10^^ h~^MQ 
in comparison to the data points in Fig. (6) of [57]. The lines show the multivariate bias (red dashed) and 
univariate bias (blue dash-dotted) non-Gaussian bispectrum as well as the Gaussian bispectrum (black solid). 
We adjust the low mass cutoff to bring the bias parameters derived from the mass function in agreement with 
the simulation measurement 6io = 2.9 and rescale &oi — > 0.75 6oi- Left panel: Bispectrum as a function of the 
high-low-fc ratio I/X3 — k/k^ and k — 0.042 /iMpc"^. Right panel: Bispectrum as a function of scale k for 
I/X3 = 7.94. 



We see that for the scales and shapes considered, the bispectrum amplitude is dominated by 
/nl and terms. The disagreement between our predictions and the simulation measurement on 
scales of fc w 0.1 hMpc~^ is probably due to loop corrections in the matter bispectrum [55] or to the 
inaccuracy of the non-Gaussian mass function that we have used to extract the bias parameters.^ 

6 Signal-to-Noise 

One goal of this work is to show the viability of the bispectrum to put constraints on primordial 
non-Gaussianity. The question is whether it is worth the additional effort of measuring the halo 
bispectrum given that the halo power spectrum can be used to put constraints on non-Gaussianity as 
well. The bispectrum analysis is naively more sensitive than the power spectrum to non-Gaussianities 
but it remains to show that this can overcome the enhanced errors. To estimate the errors we assume 
a survey of volume V from which the halo density field (5h(fc) is estimated. The bispectrum and power 
spectrum estimators are constructed from a decomposition of fc-space into spherical shells of width 
dk. Then the fc-modes within the shell are spherically averaged to obtain the estimators [59] 

7^'5(q,)^('72)<5(°ngi + 92) (6.1) 
Vv2 Jk (27r)^ Jk (27r)^ 

for the power spectrum and 

B{k,MM)^^! Pkl PvJ 7f%%i)'5(g2)5(93)^(°'(gi+g2 + 93) (6.2) 
V123 Jfci (27r)^ Jk^ (27r)^ J^^ {2ttY 

for the bispectrum, where the integrals are over spherical shells qi £ [ki — 6k/2, ki + 5k/2\. As shown 
in [59] the covariance of this bispectrum estimator is given by 

(AS)" = si23^ ('Ph(A:i) + l\ fPh(fc2) + l\ (Phiks) + (6.3) 



123 



°It is indeed possible that other non-Gaussian mass functions might give improved results, but a more careful 
treatment of these effects goes beyond the scope of our paper. 
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Figure 15. Left top panel: Variance of the halo bispectrum Bhhh for /nl ~ 100, 2:3 = 0.1. Gaussian 
(thick black) and non-Gaussian (thick blue dashed) bispectrum amplitude and error on the bispectrum for 
n = lxlO~^ft^ Mpc~^ V = 10/i"^Gpc^ kt = 2.9x10"^ /iMpc"^ and (5fc = 3fcf = 0.008 /iMpc"^ For the biased 
tracer we assume a halo sample with v = 4, corresponding to M « 1 x 10^* hT^Mq, (610 = 2.23, 620 = 0.38). 
We plot the full error Ai3(fc) (red) as well as the sampling variance (red dash-dotted) and shotnoise (red 
dashed) contributions. Left bottom panel: SNR for the halo bispectrum Bhhh (blue, red and green dashed for 
X3 = 1, 0.2 and 0.1 from bottom to top) and the halo power spectrum Phh (green dash-dotted). The signal 
is defined as the difference between the non-Gaussian and Gaussian prediction. As we go to higher k more 
squeezed configurations gain importance. Right panel: Cumulative signal-to-noise for the bispectrum (blue 
dashed) and power spectrum (green dash-dotted). We evaluate the SNR for i/ = 4 (thick) and v = 1 (thin). 



and the one for the power spectrum [59, 60] by 

(AF)' = 2^ (^Fh(fci) + (^h(fc2) + . (6.4) 

Here S123 is a symmetry factor, with S123 = 6, 2, 1 for equilateral, isosceles and general configurations, 
respectively. In the noise estimators we approximate the halo power spectrum by the leading contri- 
bution Ph(fc) = (^10 + bQi/a{k)) P\in{k) for simplicity. The fc-space volume of the fundamental cell is 
14 = fcf = {2it)^ /L^ and the norm volumes are given by 

and 

V,,= f p^j p^S^^Hq,2)^^nk,k,6k. (6.6) 

If P 3> l/n the cosmic variance dominates whereas shotnoise dominates when P <^ 1/n. 

As our fiducial case we consider a survey of y = 10 /i-^Gpc^ and 6k = 3kf = 9x10-^ hMpc^^ and 
a halo sample with overdensity 1/ — 4 corresponding to Af « 1 x 10^"* H'^Mq {bio = 2.23, 620 = 0.38). 
The number density of the halo sample is assumed to be n = 1 x 10^^ hMpc^^. The left panel of 
Fig. 15 shows the errors and the signal-to-noise ratio (SNR) for the halo bispectrum Phhh in a squeezed 
isosceles configuration k = ki = k2, X3 = 0.1. Here we define the signal as the signal-to-noise weighted 
difference between the total non-Gaussian and the fiducial Gaussian bispectrum amplitude 

SNRs(/c„/e„fc,)- (AP)2(/fc„/fc,,fcO ^^-^^ 
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with an equivalent expression for the power spectrum. We see that the /nl contribution enhances 
the signal on a wide range of scales and most relevantly on the smallest fc's. The left bottom panel of 
Fig. 15 shows the SNR for a given bin of k modes, for isosceles triangles with three different values of 
X3. We see that the bispectrum signal is enhanced as we make the triangle more and more squeezed 
and as we take the highest fc's of the triangle closer to the non-linear scale. We should stress that 
our tree-level calculation does not apply for fc's close or larger than the non-linear scale. We finally 
also plot the binned SNR from the power spectrum, which is peaked at the smallest fc's. So far 
we restricted our discussion to a particular bispectrum configuration. For the extraction of a non- 
Gaussian signal from a survey one would rather add up all the possible information, i. e. sum over all 
possible configurations up to a maximum wavcnumbcr fc,nax- The total signal-to-noise is given by a 
sum over all the possible combinations of shells up to the maximum wavenumber fcmax 



where we used fc^ = iSk. In the right panel of Fig. 15 we plot the cumulative signal-to-noise for the 
bispectrum and power spectrum for /nl = 100 and tracers with peak height 1/ = 4 and v = 1. For 
the high bias tracer we see that for the considered survey and going up to fcmax — 0.18 hMpc^^ 
we could constrain Ufj^^ 5 from the bispectrum analysis, whereas the power spectrum leads to 
constraints of ct/i^l ~ 25. Note that we are plotting the total non-Gaussian signal. To obtain the 
SNR on /nl, the lines in the figure have to be divided by /nl- For fc higher than about 0.1 hMpc~^ 
non-linear corrections should be implemented. The bispectrum wins over the power spectrum even 
if the maximum wave number is restricted to fcmax ~ 0.03 hMpc~^ . Furthermore the plot shows 
that extraction of /nl = C(10) is possible only with the bispectrum analysis but not with the power 
spectrum analysis for these tracers. Note that the choice of bias, volume and number density is similar 
to values for luminous red galaxies expected in SDSS-III (BOSS) redshift survey [61]. The bispectrum 
of the low bias v — 1 tracers shows an even more remarkable improvement in SNR compared to the 
corresponding power spectrum. 

7 Conclusions 

In this paper we present a diagrammatic prescription for the calculation of multipoint statistics 
of biased tracers of the cosmological density field accounting for the effect of non-Gaussian initial 
fluctuations. The diagrammatic approach combines biasing, non-Gaussianity and non-linear clustering 
into one consistent and intuitive picture. While we focused on the bispectrum, the generalisation to 
higher order correlators should be straightforward. Also, non-local shapes can be implemented using 
their primordial bispectrum instead of the mode coupling vertices in Fig. 6. 

We use this diagrammatic prescription to derive the halo bispectrum accounting for all tree 
level/zero loop contributions. Unsurprisingly, the bispectrum is largest in the squeezed limit 2:3 — )■ 0, 
where our prediction exceeds the results obtained using univariate bias only, by about a factor of 2 on 
scales of fc ss 0.03 hMpc~^. Given the fact, that we see these corrections compared to the univariate 
bias approach already at tree level, we caution the use of the univariate bias expansion as presented 
in [22] and [21]. While the qualitative results are quite general, the quantitative results presented in 
Section 5.3 are somewhat dependent on the choice of the halo sample. We considered a i' — A tracer, 
corresponding to haloes of mass A/ sa 1 x 10^^ H^^Mq {bio — 2.23, 620 — 0.38). 

Comparing our results to the ones obtained in [21] we see that their one-loop terms proportional 
to 62 are comparable to some of our tree level terms if i) the smoothing scale equals the halo mass 
scale, ii) high mass haloes are considered and iii) large scales are considered. But, these terms can 
not be seen as a replacement to the non-Gaussian terms obtained in our approach because all the 
one loop terms arising in the univariate bias approach are also present in our approach. The final 
comparison of the terms would require a careful resummation of the bias parameters, which was not 
considered in [21]. 




(6.8) 
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Probably the most important result of this paper is that the halo bispectrum analysis offers an 
alternative to the power spectrum for detecting local non-Gaussianities with an even higher constrain- 
ing power. For the i> ^ 4 tracers in a = 10 ft.^'^Gpc'^ survey, our signal-to-noise analysis predicts 
a factor of five improvement in the constraints on /nl compared to the power spectrum. For lower 
bias tracers at the same number density the total signal-to-noise is a bit lower, but it is in fact much 
higher relative to the power spectrum analysis, which contains no signal for unbiased tracers. 

Our results suggest that the bispectrum should be the statistic of choice for detecting primordial 
non-Gaussianity in current and future survey data. However, some additional work has to be done 
before this method can be applied to the real data. One extension is to make predictions for pho- 
tometric surveys, where only projected density fluctuations are observed. Alternatively, if a rcdshift 
survey is used in the analysis then the predictions here should be generalized to include redshift space 
distortions [62] . In the light of ever increasing surveys and simulations one might also be concerned 
about general relativistic corrections on horizon scales [33] . 

Convergence is probably one of the most important problems for perturbative calculations. For 
Gaussian initial conditions, comparisons of the matter bispectrum in simulations to the theoretical 
predictions as presented in [55] and [63] conclude that 1-loop calculations are required to achieve 
a reasonable agreement on scales of fc « 0.1 h~^Mrpc. The bispectrum in presence of local non- 
Gaussianity is the perfect statistic to apply perturbation theory at tree level combined with the bias 
from the peak-background split since the non-Gaussian effects are most prominent for low k or large 
smoothing scales. Still, we focus on a tree level calculation and find that the signal receives relevant 
contributions from scales close to the non-linear scale. For this reason, loop-corrections should be 
examined to fully assess the detailed amplitude of the signal. As we discuss, this will probably require 
to study the cutoff dependence both of the loop corrections and of the bias coefficients, such that 
final observables do not depend on the cutoff.^ Some of these 1-loop terms will be renormalized and 
absorbed into the tree level terms discussed here. Still, a consistent calculation of the halo bispectrum 
at the next order would require consideration of non-linear couplings up to F4 and biasing up to fourth 
order, both in 6 and ip. This increases the number of terms by a large amount and goes beyond the 
scope of the current paper. The final decision about the validity of the perturbative calculation 
presented in this paper has to be based on a detailed comparison to the halo bispectrum measured 
in simulations. A first comparison of our results to the measurements published by [57] shows an 
encouraging level of agreement. 
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A Standard Perturbation Theory 

This appendix reviews the essence of cosmological perturbation theory and serves as a source for the 
most important equations. For a more detailed treatment we refer the reader to the comprehensive 
review on the subject by [47]. The evolution equations for the cosmic fluid in an expanding Friedmann- 
Robertson- Walker Universe can be formulated in terms of the overdensity 5 and the velocity divergence 
6* = V • t> as 



d9{x,T) 
dr 




(A.l) 



(A.2) 
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A direct solution to this coupled differential equations does not exist. Therefore one has to restrain 
to a perturbative solution, expanding the density field and velocity divergence in a power series 



^?(fc) = -H^0(")(fe), 



(A.3) 



where S^\> and 0^"'' are 0{S^ p). The solutions for the n-th order contribution to the fields are given 

by 

d3 



(27r)3 



e^''\k) = J^---J ^W(gi)---^m,p(q„)G„(qi,...,q„)'5'°^('7i + --- + gn-fc), (A.5) 

where the coupling kernels Fn and G„ can be obtained using recursion relations. For the results 
presented in this paper we need Fi = Gi = I, and 
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Using the above results, the one-loop corrections to the matter power spectrum read as 



P22 = 2 



(2^)^ 



:P{q)P{\k-q\)\F2{q,k-q)\\ 



and 



PlZ = QP{k) J ^P{q)F,ik,q,-q). 



(A.6) 

(A.7) 
(A.8) 



B Explicit Derivation of the Bias Parameters 

As shown in Section 3, long wavelength modes in presence of primordial non-Gaussianity change the 
effective collapse threshold and variance of density fiuctuations. In the mass function, the presence 
of a long wavelength mode can be accounted for by rescaling the peak height as 



Sc - Si 

ctg(1 + 2/nl<^i + igNLff + 



(B.l) 



As we will need them for the explicit calculation of the bias parameters, we write down the par- 
tial derivatives of the peak height with respect to primordial potential and long wavelength density 
fiuct nation 
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ddf 
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/NL- 



5^ 



(24/2^ - 12gNL)z^ 



(B.2) 
(B.3) 

(B.4) 



Using these expressions and the results of Section 3 we can first derive the Gaussian bias parameters 
arising from the rescaling of the collapse threshold 



1 dn 1 dn dv 



1 2v dn 
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0. 



^° ndSi^ ndv^ \d6j n dv dSy^ 

4 u"^ d'^n ^2 u dn (B 6) 

n (5^ dv"^ n dv 

The non-Gaussian bias parameters are now arising from the rescahng of the variance in the denomi- 
nator of Eq. (B.l) 

Oqi —--^ — --^-^ — — = ^ - ^/nlOcOio [tij) 

n oipi n av Oip\ n ov 

1 d^n 1 d^n dv dv 1 dn d^v 
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where we used the results for b^Q and b^Q obtained above to replace the partial derivatives of the 
massfunction with respect to the peak height. 
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